from math import sin,cos,pi,sqrt,e

def gabor_function(x,y,envelope_rotation,sinusoid_frequency,sinusoid_rotation,envelope_height,envelope_width,envelope_peak_x,envelope_peak_y,sinusoid_offset):
    def complex_sinusoid(F,omega,P,x,y):
        u = F*cos(omega)
        v = F*sin(omega)
        return cos(2.0*pi*(u*x+v*y)+P)
    def gaussian_envelope(x0,y0,a,b,phi,x,y):
        x,y = xphi(x,y,x0,y0,phi), yphi(x,y,x0,y0,phi)
        return e**(-pi*((a**2*x**2)+(b**2*y**2)))        
    def xphi(x, y, x0, y0, phi):
        return (x-x0)*cos(phi)+(y-y0)*sin(phi)
    def yphi(x, y, x0, y0, phi):
        return (y-y0)*cos(phi)-(x-x0)*sin(phi)    
    return (
        gaussian_envelope(
            envelope_peak_x,
            envelope_peak_y,
            1.0/envelope_width,
            1.0/envelope_height,
            envelope_rotation,
            x,
            y
        )
        *
        complex_sinusoid(
            sinusoid_frequency,
            sinusoid_rotation,
            sinusoid_offset,
            x,
            y
        )
    )
    
def gabor_filter(width,height, envelope_rotation = 0.0, sinusoid_frequency = 1.0/20.0, sinusoid_rotation = 0.0,envelope_height = 80.0,envelope_width = 80.0,envelope_peak_x = 40.0,envelope_peak_y = 40.0,sinusoid_offset = 0.0):
    return [[gabor_function(
            float(i),
            float(j),
            envelope_rotation,
            sinusoid_frequency, 
            sinusoid_rotation,
            envelope_height,
            envelope_width,
            envelope_peak_x,
            envelope_peak_y,
            sinusoid_offset
        ) for i in xrange(width)] for j in xrange(height)]